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In the current article, we rederive the lattice-fluid excess models UNIQUAC, UNIFAC, and COSMO-RS 
from a continuum functional. The calculation explains the missing dependence on the particle geometry and 
how to include the Coulomb interaction, problems that are common to all three models. It is then shown that 
the Wilson ansatz, used in UNIQUAC and UNIFAC to minimize the grand potential, is not a physically valid 
solution of the Euler-Lagrange equation. A consistent approach is the Larsen-Rasmussen equation, which 
forms the foundation of COSMO-RS. We then analyze the various approximation methods and interpret 
them in the framework of a molecular density functional. 
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1, Introduction 

A fundamental problem of density functional theory (DFT) is the construction of the 
grand-canonical potential, which in soft-matter physics has to cover the wide spectrum of 
interparticle potentials ranging from the strongly repulsive inner domains of molecules to 
the Coulomb interaction of ions HI IS . A common strategy is therefore the separation of 
scales and to approximate the repulsive inner core by a hard-particle potential. Assuming 
that the residual interactions are weak, this reduces the construction of the grand potential 
to the perturbative expansion of the free energy and the derivation of the distribution 
functionals for hard particles. The latter have been derived recently, which leaves us to 
identify a suitable representation for the molecular functional 

For practical applications, it is necessary that the minimization of the grand potential 
remains calculational manageable. We therefore require the functional to be of low pertur¬ 
bation order, valid for short and long-range interactions, to allow further approximations 
in the hard-particle functionals, and to be numerically simple to integrate. It is not easy 
to meet all these criteria. But a good starting point is the analysis of lattice-fluid excess 
theories, especially the well established universal quasi-chemical model (UNIQUAC), its 
extension to functional activity coefficients (UNIFAC), and the conductor-like screening 
model for real solvents (COSMO-RS), which are simple but fully operational density 
functionals that derive the liquid-liquid equilibrium by minimizing the grand potential 

M- 

Lattice theories were among the earliest models in solid-state physics, best known for 
the Ising model and its solution in one and two dimensions |[9l-[TTl . But they were also 
applied to the liquid state by Flory, Huggins, Staverman, Guggenheim, et al., until ex¬ 
perimental results proved the more random structure of the particle ordering l[T2] - [T5l . 
As a result, they lost its popularity in physics, but remained in use in biology, chem- 
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istry, and engineering, where they have been extended hy Prausnitz, Ahrams, and Mau¬ 
rer to the UNIQUAC, UNIFAC models and hy Klamt to COSMO-RS, with the alterna¬ 
tive implementation of COSMO-SAC hy Lin and Sandler ll^ [T6] - l22l . Together with the 
group-contrihution approximation, they successfully predict the liquid-liquid equilibria 
for a large class of molecules using a simplified interaction potential for their chemical 
groups or surface charges. 

A major disadvantage is their limitation to a fixed value of densify and pressure as a 
consequence of fhe laffice represenfafion of fhe molecules by linear chains of laffice cells. 
This allows fo disfinguish individual parficles by fheir posifions, which resulfs in a wrong 
combinaforial facfor of fhe free energy 0. The laffice functional ifself is fherefore ill 
defined. Buf fhe incorrecf parf cancels in fhe excess energy of a mixfure wifh parficles 
of maximal packing fraction. This parfly correcfs fhe error of fhe laffice ansafz, buf also 
shows fhaf any improvemenf of fhe model is only possible in fhe framework of DFT. 

Sec. [^begins wifh a general discussion of fhe free-energy represenfafions and fheir re- 
specfive perfurbafive expansions. If is fhen shown in Sec.j^fhaf fhe laffice models originafe 
from fhe dual functional. We analyze fheir solufions of fhe Euler-Lagrange equation and 
reinferpref fhe underlying approximations in ferms of a continuum molecular functional. 


2. The dual Free-Energy Functional 


The Hohenberg, Kohn, Mermin fheorem proves fhe unique mufual relafionship befween 
interacfion and grand-canonical pofential IB, forming the foundation of DFT, as it guar¬ 
antees that different representations of the functional determine the same thermodynamic 
ground state. It thus constraints the number of alternative representations for a given in¬ 
teraction potential, which have to be related by similarity transformations. Ignoring map¬ 
pings that correspond to internal symmetries of the potential or result in contributions 
which cancel in the Euler-Lagrange equations, the only nontrivial symmetry is the Leg¬ 
endre transformation of the grand potential D., exchanging its canonically conjugate vari¬ 
ables. Eor the simplest case of pair interactions, it replaces the potential ^ij by its dual 
pair-correlation functional gij, defined by 


do. 1 

d^j ~ ■ 


( 1 ) 


This shows fhat ^{<pij) has only two representations, either as the free-energy functional 
Q.^{(j)ij) or its Legendre-dual counterpart QP{gij). 

Most molecular functionals use QF as the starting point, as its perturbation expansion 
in r-particle densities Piy...i^ is algebraically well understood O]. But it will be shown in 
the next section that the lattice models derive from the dual functional whose ana¬ 
lytic form is more complex and the perturbation expansion of g 2 does not result in either 
the direct or the distribution functionals. Lor comparison, we will derive both representa¬ 
tions, analyze their respective perturbation expansions, and discuss their advantages and 
limitations. 

Beginning with the free-energy representation, the functional (/),y) of 

the particle density pi depends on the inverse temperature jS = 1 /k^T, chemical potential 
Pi, and external potential (pf^ for a mixture of / = 1,... ,M compounds. It is an integral 



[p/(ln (pA) - 1) - PpiiPi - PH] dji - coipi) 


( 2 ) 


over the positions and orientations y;- G M” x SO(n) of the «-dimensional Euclidean space 
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and depending on the thermal wavelength A, and direct correlation functional co(p/). 

Its perturbation expansion in the potential (p = (p^ (p^ of hard-particle <p^ and soft 

interaction (p^ is a formal Taylor series of the logarithm of the grand canonical partition 
integral, whose first and second order in the Mayer functions /I have the form 


_ 1 /"foH. . . H ws .S Jy. . . . 


(3) 


where a sum over paired indices is implied ||T1 . Higher order terms are readily obtained 
using diagrammatic techniques Q, where each product [/f ]”* couples to a homogeneous 
polynomial of r-particle densities of order 2 <r < 2m, integrated over the {m — l)n{n + 
l)/2 coordinates of position and orientation. The rapid increase in the dimensionality of 
the integrals effectively limits the perturbative approach to the first or second order. 

The expansion ([^ requires the particles to interact by pair potentials. But the same 
approach also applies to irreducible m-particle interactions, when the fully /|-bonded 
subdiagrams [/f ]'” are replaced by the Mayer function For the 3-particle interaction 
^,i,' 2 / 3 , e.g., this adds the leading correction 


PDF = PDF{(Pi^ 


5 J Piihh-Fiihh ^y^hh + • 


(4) 


to the 2-particle functional QF{(pij). 

The second representation is the dual functional First derived by Morita and Hi- 
roike using diagrammatic techniques ll23l - lT7l . it replaces (pij by its canonically conjugate 
variable gij. To perform the Legendre transformation, we integrate ([^ over d(pij 

n = + -Jp.PjgijfpijdYij - - JPiPjPij5g,jdYu . (5) 

To complete the integration over 5gij, Morita and Hiroike derive the self-consistent clo¬ 
sure relation between (pjj and gij lU : 

In(gp) = -Ppij+dij + hij-dj , (6) 

introducing the bridge functional dij of 2-path connected clusters, the pair correlation 
hij = gij — 1, and the 2-particle direct correlation functionals cij. To eliminate the re¬ 
maining dependence on the free-energy representation, aj is then substituted using the 
Ornstein-Zernike equation 

oo 

= I 

n—Y 

Inserted into Q, they form the infinite sum over h-bonded ring integrals 

“ r (—1)” 

PiPj{Cij ~ hij) ShijdYiJ = / Pki ■ • •Pk„hkik2 - • •^k„kidYki...k„ (8) 

n=3“' ” 

in the final representation of the dual grand-canonical potential |[2^ \2M : 

= £ / pi\n{piKi) - Pi - ppiPidYi + \Y. f P‘Pj{siA^iSij)-8ij + ^)dYij 

i ^ ij 


j { Pki ■ ■ ■ Pkn ^iki ■ ■ ■ hk„J dYki...k„ ■ 


(7) 
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B f 1 “ /■ (—1)” 

+ 57 52 / PiPj Sij *Pij dji] + 7 52 / ~ Ph ■■ - Pkn ^kik2 ■ ■ ■ ^k^ki djki ...k„ (9) 

^ ij ■' ^ n=3-' " 

PiPjdij^SijdYu. 

where the integration constant +1 in the second integral has been chosen to reproduce the 
ideal gas in the limit (pij —i- 0. 

Compared to the free energy representation Q, the analytic structure of the dual func¬ 
tional is considerably more complex, while containing exactly the same information for 
pairwise interacting particles. A common simplification is therefore to set dij = 0 and to 
use either the Percus-Yevick (PY) or the hypernetted chain approximation (HNC) for the 
closure relation Q 

PY : gij exp (j3 j) = exp {dij + hi - c/y) « 1 + hij - aj 
HNC : In (gij) ^ -I5(l)ij + hij - aj . 

In combination with the Ornstein-Zernike equation HI, they provide easier to solve self- 
consistent integral equations for h 2 and C 2 - 

Probably the best known example is the PY approximation for hard spheres and its 
solution for g 2 developed by Wertheim, Thiele, and Baxter ll28l - l3Tl . Another example is 
the Coulomb potential (j) = q^jr for point-particles of charge ±q. Its slow radial decline 
allows the long-range approximation C 2 for which the HNC equation can be 

solved, using the Fourier transformation C 2 = F{ 02 ) to decouple the Ornstein-Zernike 
equation 

/V 2 

ln{g 2 ) = -l5^ + h 2 -C 2 ^h 2 =F^^( ) =-j3 —exp(-k£,r) . (11) 

VI — pc 2 / r 


This result reproduces the characteristic Debye-Hiickel screening for an ionic liquid of 
wavenumber ko = {ATi^pq^y/'^ and, together with the infinite sum over the ring integrals 

JP (q - h2)5h2 df = ph 2 - ^P^hj - In (1 + pk) , (12) 

yields the Debye-Hiickel functional for charged spherical particles 111 |32l . 

The calculation illustrates how the combination of Ornstein-Zernike and closure equa¬ 
tion Q improves the low order approximation. Actually, it is an example of a duality 
transformation that inverts the length scales by exchanging a pair of canonically conju¬ 
gate variables, mapping the perturbative sector of one functional to the non-perturbative 
of its dual. This shows that the two representations and although equivalent in 
their total information, have different application ranges when the perturbation series are 
restricted to a finite order. 

Contrary to the expansion O, the perturbation theory for is an expansion of g 2 in the 

soft correction term Fij = of the Mayer fij = + Fij and Boltzmann functions. 

Its lowest order diagrams are shown in Fig. [T] illustrating the successive replacement of 
-bonds by f-functions. Its corresponding substitution in the functional is generated by 
the formal derivative 


Sij = elgfj + el 


8g^j 

SfH 


Fikdjk + 2^0 


~ gi}\2fl+ gii\2.\ FgijUp + • ■ • , 


TYH ^hh dTkik2 + • • • 
^Jkik2 


(13) 
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a) 


b) 


S2 = ° .° 




c) 




+ 



Figure 1 . The perturbative expansion of the pair-correlation functional g2 written in Mayer diagrams: solid lines indicate 
/^-bonds, dashed lines correspond to F-bonds, and dotted ones represent the product of hard- and soft- Boltzmann 
functions. The main contribution derives from the fully bonded diagrams a) with no internal F-bond g2\2.o, b) one F-bond 
linked to a rooted node g212.1. and c) one internal F-bond g2 14.1. 


whose contributions can be further approximated by a series of correlations g2\n,k of 
hard-particle intersection centers and k internal F-bonds. Using the notation of H |3, 
the functionals of lowest intersection order and with at most one internal F-bond have 
intersection diagrams of the form: 

• e- ■ +e- ■ 

'''■2.1 • lb 

^' 3,1 • ^ni2^i2i3 la lb '' '' 

^' 4,1 • ^'1'2'^0'4/a lb Ic Id ’ 

whose resummation yields the functionals 


h h 

gilhhfl = 4ii2^hi2i^AI + %^b . J" ) 

A ^ab 


(15) 


with the correlation integral of the inner bond 

yabcd = jFi3iAi^'>'cPi3)i^'b^'dPiA)d'yi3H ' (^ 6 ) 

This example shows that the representation in intersection centers can also be extended 
to the perturbative expansion of g 2 , completing previous results for the direct and dis¬ 
tribution functionals SO. But for many applications it will be sufficient to restrict the 
series to the first order in the hard-particle correlations and to use g 2 ~ ^ 2^2 1 2 , 0 - This 
approximation has the additional advantage to satisfy g 2 ^) > 0 and g 2 {r ^ °°) = ^ for 
any r G M'’. These two conditions are easily tested and provide an important advantage 
over where no such constraints exists. Because — 1 < /f < °° and the discontinuity of 
§2 at particle contact, the integral of its first approximation order ffg^ is indefinite and 
very sensitive to changes of the hard-particle geometry and soft-interaction potential. An¬ 
other disadvantage is the dependence of ^ on Cq and g^, which requires the redundant 
calculation of two functionals with the same information content. 

In summary, depends only on one class of correlation functionals, allows a sim¬ 
ple consistency test for its perturbative corrections, and applies for short- and long-range 
interactions alike. A disadvantage, however, is its limitation to pair interactions. But as 
higher-order irreducible m-particle potentials are often very short-ranged, they can be 
coupled perturbatively using the expansion Q. Thus, despite its complex structure, the 
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dual functional is the preferred ansatz for the perturbative construction of a molecular 
model, which will he further supported in the next section, where we make contact with 
the lattice theories. 


3. Lattice-Fluid Models derived from the Dual Functional 

Lattice models for fluids use a discretized Euclidean space, with molecules represented hy 
linear chains of cells, reducing the configuration integral to a sum over all allowed particle 
insertions |l7l|8j[TT]|T5l. To simplify the derivation, two additional assumptions are made: 
1. molecules are closely stacked, i.e. the packing fraction for all systems is T] = 1, and 2. 
interactions only occur between next neighbors. 

A mixture of N = particles with i = 1,... ,M compounds is therefore independent 
of volume effects, from which follows that the excess free energy of mixing 


F^{{xi}) = F{{xi}) -Y,XiFi{xi = 1) 


(17) 


is a function of temperature and the molar fractions Xi = Ni/N. 

To derive the lattice model from the functional we have to interpret its two assump¬ 
tions in terms of the continuum formulation. The first constraint of close packing is readily 
implemented for a mixture of constant densities p, = A, /V and their pure-compound sys¬ 
tems Pi = Ni/Vi with partial volumes V) = xiV and molecular volumes v,-: 


I = rj = Y^PkVk = PiVi = f)i . 


(18) 


k 


The second constraint reduces the correlation length of the pair-distribution function to its 
next neighbors. In a first step, we therefore neglect all g 2 in the functional Q beyond the 
leading order 



i •' ij 


thus removing the bridge and ring integrals responsible for the Debye-Hiickel screening. 
Next, we restrict the spacial range of g 2 to the first particle shell A,j, which comes closest 
to the idea of next-neighbor correlations between cell elements. Introducing the definitions 


j ‘j ‘j 

the pair density and correlation functions can be rewritten 




( 21 ) 


Zi = 



in terms of the local coordination number 6/y = Zijjzj, surface fraction 6/ = PiZijz, and 
the volume fraction of lattice theories. These variables are not independent, but related 
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by the permutation symmetry pij = pji and normalization 

dijdj = djidi, Y. f 1 ’ L / 1 • ( 22 ) 

i •'^ij i •'^ij 

Using these identities, the continuum functional ( [T^ can now be written in the basis of 
lattice variables. 

Beginning with the potential energy 


/ PiPjSij^ijdjij = I PiPjSij^ijdjij = ^zNYPijPjP 


IVxV 


ij > 


(23) 


the integration over V xV separates into a sum over N/2 particle pairs of volume A/y, 
while the potential (pij is approximated by the constant energy parameter Eij of neigh¬ 
boring cells molecules. The same transformation also applies to the logarithmic term of 


tE ( PiPj8iM8ij)dYij = 

z JVxV z ■■ \ tii z ' 

(24) 

= lz«p,;%ln(|')+z«£e,ln(|) + izWl„(^), 

whose constant contribution cancels in the excess free energy ( [TT] ). The same applies to 
the linear term 


1 

2 



PiP}{8ij-^)dyij = 


^A(2z-A). 


(25) 


Slightly more complicated is the transformation of the kinetic energy, as the integration 
over the configuration space A,j effectively reduces the number of independently moving 
molecules. The excess kinetic energy of unpaired particles 

= E / (pi 1*1 (Pi^i^ -Pi)dri-Y f (pi 1** (A'^') “ A) dji 

i '^V I J Vi 

(26) 

= EA^/ln(^) =aEjc;I n ((/>/) , 

has to be corrected by the kinetic energy of the particle pairs. To determine their contri¬ 
bution, one has to observe that the translational and rotational degrees of freedom of one 
particle is bound to the second particle of the cluster. We therefore subtract for each pair 
the kinetic energy of one particle. The amount of energy bound by the density of p,z//2 
particle pairs is determined by the difference: 


\'Ljy(Pi^i^^(Pi^i^i)^Pi^i^dri-^Yjy^i(Pi^'^(Pi^i^^Pi^dri 


( 27 ) 
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From this derives the excess kinetic energy stored hy the clusters 


= 0 L / Pi^i Pi^i 


Oi z 




Pi\ 




Subtracting this result from the excess kinetic energy of the free particles 
effective kinetic energy of the lattice fluid 


( 28 ) 


yields the 


1 


In m - xzA^E In 


(29) 


Combining the previous results (231, (251, (24 1 with the identities da = 1, 6, = 1 for the 
pure compounds and canceling constant contributions, we finally arrive at the excess free 
energy of the lattice liquid 


PF^/N = In (,/>,) + I ^ ean ( 2) + IE ( ^) + ^ (% - 


6i\ ^ _ r. /ft 


ft 


(30) 


whose first two terms correspond to the results from Flory-Huggins and Staverman- 
Guggenheim |[T2] - [T5]| . The corresponding grand-canonical excess functional follows by 
adding the excess chemical potential of paired particles 

a^idp) = F^{e,j) - zivE ^iPi ■ (31) 

i 

The mixtures are now uniquely determined by the four sets of variables dij, ft, ft, z, 
and the constraints (22 1 . But only dij is fixed by the Euler-Lagrange equation of 
The remaining variables still need to be determined from their definitions (20), (21), and 
the assumptions of the lattice model. The molecules, e.g., are flexible, linear chains of 
cells without self-intersection. Their specific shape is therefore undefined, but the volumes 
Vi and surfaces a, are constant and the contact probability independent of positions and 


orientations, thus corresponding to g 
function 


O'I A,7 


i c. The hard-particle pair-correlation is then a 


sfj{t)\Kij = ceij{t)5{t) 


(32) 


for particles whose surfaces are separated by the distance t = 0, as shown in Fig.[^). 

To calculate the coordination numbers (20), we use the representation from App.[A]for 
the integral measure dji] of two particles with principal curvatures at a distance t = 0 
and rotation angle 0 < (p < In. Expanding the determinant (A 10) 

det[A^^^ + -I- -|-sin^ (^)(fc{^^fcp^ -f 

-fcos^ -I- 


and integrating (32) over all relative orientations of the two particles, yields the surface of 
the Weyl tube: 



Stt^c J efj{t)8{t)dtdOi = %n^caidij 


(34) 
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and the surface of their Minkowski sum: 


sfjdjij 


:J det[A*'^^ + M ^X^^^u]d(l)daidaj = ?iK^c{ai + aj + . 


(35) 


The Minkowski surface is therefore not simply the sum of its individual surfaces hut cor¬ 
rected hy the product of mean curvatures . Its counterpart in the lattice representation 
are cell segments adjoined at the edges of the molecule hut not on its surface segments. 
These cells, however, are ignored in the next-neighhor approximation, explaining why 
the lattice models cannot represent the specific geometry of a particle. Omitting the non¬ 
additive part and introducing the packing fraction rj = PiVi, determines the remaining 
three groups of variables 




XjVj 

LkXkVk ’ 


ei = 


XjUj 

^kXk^k 



(36) 


as a function of the universal model parameter zq- 
The last, hut subtle, step in determining the thermodynamic equilibrium is the mini¬ 
mization of the functional 


5QP 


5QP iSOPdgij^^ ^ 

bPk 82 2 bgij bpk 


(37) 


The Euler-Lagrange equation of the first term defines the chemical potential, while the 
second reproduces the constraint Q as a self-consistent equation. To compare this equa¬ 
tion to its analogue in we apply the previous approximations by omitting terms of g 2 
beyond the linear order In (§ 2 ) = —jS <p + d 2 + h 2 — C 2 ^ — j3 (j) and rewrite the correlation 
function in the basis of the lattice variables (|2T]) 


SijlXij — 


^ijZiZj 

di z 


exp {-l5(j)ij)\Aij = exp i-pEij) 


(38) 


The corresponding minimization of with respect to dij and the constraints (22 1 yield 
the Euler-Lagrange equation for the lattice model lITSl 


5^^ 


= 0 : 


Ondjj 


— exp ( l5^2,Eij Eii 


(39) 


for which we introduce the notations: 


T,-,:—exp( P[2eij Eu £p]) , by:—exp ( £ 77 ]), '^ij — Ujtji 


(40) 


By inserting (38 1 into gijgji/{giigjj), it is easily seen that the minimum of the continuum 
functional and that of its lattice counterpart ( [39l ) agree, therefore proving that the first-shell 
approximation does not violate the thermodynamic consistency of the functional. 

In the literature, two different solutions for ( |39l ) can be found. The first one, developed 
by Larsen and Rasmussen (LR) |[33l . uses the symmetry properties (22 1 to derive the 
square root 


(Pij\^ Pa Pjj 2 12 Pa 


1 

I’j 




(41) 
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which can he numerically solved for bi and hack-inserted to obtain Bij. The alternative 
approach goes hack to Wilson (W) ll^ and uses the ad hoc separation 


a ~ p ‘‘j ^ Q ~ a 
tijj Uj Uii «/ 


d.j = 


^ij 

Hk tk j 


dj, = 


^j^ji 
Hit ^k tki 


(42) 


to obtain two independent solutions for (391 in terms of tij. This approach, however, is 
inconsistent, as can be seen from the missing permutation invariance of tij in its indices 
and by inserting (38) into gijjgjf. 


- tif 

Ojj Oj zi • 


(43) 


The Wilson ansatz is therefore only a formal solution, depending either on the volume or 
the surface fraction and at most applicable for molecules of similar spherical size Zi ~ Zj- 
Inserting ( [44] ), ( |4^ into the functional ( |30| ) and taking account of the two independent 
solutions of the Wilson model, yields the minimum of the excess free-energy with respect 
to dij 


In (</>/) + IL In (e,-In 


L ft J ’ 


PFl/N = f Eftln ( J) -zA^Eftln 

i i ^ ri i i 


(44) 

(45) 


where the second result corresponds to the UNIQUAC model introduced by Prausnitz, 
Abrams, and Maurer ifT^FTTll . 

The liquid-liquid equilibrium at a given reference point of density and pressure is now 
determined by the excess free-energy function and the parameters v;, a,-, T/y and t,y re¬ 
spectively. Their values can be adjusted to experimental data if a sufficiently large set 
is known. This is especially convenient for the analytical solution of the Wilson ansatz 


(42 1 , which partly explains the popularity of the UNIQUAC model. If, however, the data 


set is too small, one has to resort to further models to specify the geometry and inter- 
molecular potentials. One such approach is the group-contribution approximation, which 
uses the observation that the chemical and physical properties of organic compounds are 
often dominated by their functional groups. Together with the lattice assumption of next- 
neighbor interactions, the potential fty is replaced by a superposition of interactions ^^^3 
of its a, j3 = 1 ,... ,Nc functional groups 


<^ij =L 

aj 3 


klj (f>al3 ) 


(46) 


related to an analogous transformation of the pair-correlation functionals 

5a 

afi 


5(^i 


5a 




Pi] — ^PiPjkia^pSap ■ 

ap 


(47) 


The functional groups are the lattice equivalent of the site-site interactions used for 
molecular fluids lIH. But in combination with the next-neighbor approach, they decou¬ 
ple and formally replace the molecules as individual particles in the potential part of 
the free energy. Writing its contribution in group indices, the transformation leaves the 
particle density p,- and the product of canonically conjugate variables gqfty = gap'Pap 
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invariant. Only the integral measure djij = nfn^jdjafi is changed hy the transition dOi = 
\dOi/dqa\dqa from surface elements to surface groups or charges qa- The transformation 
of the potential energy thus remains formally invariant 

E / PiPjSij<i>ijdYij = Y.Y. / PiPjSafi<i>apnfn^jdYa^ = E / PaPpgaP^apdYap , (48) 

ij ij aji-^ a/3 


if the particle density is redefined as the density of group elements 

Pa = • 


(49) 


Because the partial integration Q commutes with the coordinate change (46 1 , the same 
transformation applies to the complete functional. The excess free energy and the 


models (j4^, ( [45| ) therefore remain formally invariant. Using the substitution x, = 
for the molar fractions, yields the lattice variables of the group-contribution models 


n?Xa 


Y^iXinfag 

Yk^k^k 


Paj8 ) ^a/3 j fa/3 j 


(50) 


for the group surface ag = n‘^ai and the group volume Vg = n'gVi. 

Writing the UNIQUAC equation in the basis of group contributions reproduces the 
UNIFAC model llTSl . Its extended class of parametrized functional groups improves the 
accuracy of the UNIQUAC model and allows to interpolate between molecules of similar 
chemical structure. But its dependence on the Wilson ansatz, the low spacial resolution of 
the interaction potential, and the heuristic notion of functional groups limits its value as a 
guideline for further improvements. 

An approach that avoids these complications is the COMOS-RS model |0|T3. Instead 
of the functional groups it uses partial charges qi^ localized at the segments of the 
discretized surface of the molecule. Their values are derived by a quantum mechanical 
COSMO calculation, approximating the dielectric background of the liquid by the bound¬ 
ary condition of a conducting surface, which can be solved by inserting mirror charges 
—lESl. Removing the boundary condition, these charges generate an electrical field E 
outside the particle, pointing into the normal direction of the segments . This has 
to be taken into account, when the interaction energy between two neighboring molecules 
is determined. Using the Maxwell tensor Oab = 1 /(47r) — E'^dabl'Yj, the energy of 

the electric field between the surface charges qi^, qj^, separated by a distance is 
approximately 




Knjp 

^‘ajp 


(51) 


Inserting this result into (40), yields the interaction matrix for ([44]) 


'-tajp 


= exp [- y — [qia nia - qjp njp) ] 


'ajp 


(52) 


Solving the self-consistent equation is still a time-consuming task even for small 
molecules. Given a mixture of particles with S/ surface segments, the rank of the ma¬ 
trix is Si -|- ... -I- Sm, which for binary mixtures is of order 10^ — 10^. To shorten the 
calculation time, the COSMO-RS model introduces group variables to coarse grain the 
number of charges qg = n‘^qi^ and segments ag = lYgai^, simplifying the self-consistent 
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Figure 2 . Comparing the interaction models of COSMO-RS and the molecular functional: a) The COSMO-RS model 
maps the surfaces of the molecules to unit spheres, with partial charges interacting over a fixed distance tQ. b) The molecular 
functional determines the pair-correlation functional, coupling the hard-particle geometry to the soft interaction of partial 
charges. The grand potential is the integral over all segment-segment combinations, distances t, and axial rotations (j). 


equation |@ 


^ = E Oaba , = exp [- ^ K{qa + )^] , (53) 

for molecules of antiparallel surface segments ni^ = —fij^ and separated by an average 
distance = to, whose value has been absorbed in the overall constant K. Together 
with the reference geometry of the unit sphere, this corresponds to the interaction model 
shown in Fig.|^). 

Apart from the electrostatic interaction, the COSMO-RS model also includes disper¬ 
sion effects and hydrogen bonding, but fails for the Coulomb interaction. This is to be 
expected, as the next-neighbor ansatz requires the correlation length to be of the order of 
the first particle shell, whereas the correlation length of strong electrolytes is significantly 
larger. As a result, the ring integrals Q can no longer be ignored. A first approximation is 
therefore to couple the Debye-Hiickel ( [T^ or Debye-Hiickel-Pitzer term to the grand po¬ 
tential and to derive the new closure condition from the Euler-Lagrange equation. 

For the example of an electrolyte with point charges ±q, this yields an implicit equation 
in g 2 


ln(g2) = -P(l> -F 


A. 

Phj \ 

l+phi' 


(54) 


Its algebraic solution is no longer possible. But the pair correlation is still dominated by 
the potential ^ at particle contact r = to and only modified by fhe Debye-Hiickel ferm 
af disfances ro = iTz/ko- If ro ~ to, the screening effect is small and can therefore be 
ignored. Whereas the detailed geometry of strong electrolytes ro »to is less important for 


larger distances, and the pair-correlation function /i 2 (ro) can be approximated by (11 1 . The 
inhomogeneous charge distribution has then the effect of a charged dielectric background 
that contributes a compensating potential to ^ without violating the additive structure 
required for the self-consistent equation of the lattice model ( [4T] ). 

Despite this generalization, lattice models remain limited by the fixed reference value 
of density and pressure and the neglect of the molecular geometry. Any improvement re¬ 
quires the construction of a density functional. A natural link between both descriptions is 
the COSMO model. The cavity and partial charges provide the necessary information to 
define fhe hard-parficle geomefry and soff inferacfion for fhe approximafe pair-correlafion 
funcfional gij{Oi,Oj,tij,<^), as shown in Fig. [^). If is a funcfion of fhe surfaces d/, Oj, 
separafed by fhe disfance tij, and rofafed by fhe axial angle (p. If also infroduces corre- 
lafions befween spafially separafed surface segmenfs, which is unavoidable fo describe 
elongated or concave molecules and fo approach problems from biology and chemisfry. 
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4. Discussion and Conclusion 

The last two sections have shown that is a promising starting point for the construction 
of molecular density functionals. It has a simple first perturbation order, a local consis¬ 
tency test for the pair-correlation functional, and it reproduces the excess free energy of the 
lattice models. By combining the various approximation methods, we have now a better 
understanding for the continuum functional and how to combine the quantum mechanical 
data of a COSMO calculation with the simple interaction model of classical mechanics. 

On the other hand, the derivation also gives new insight into the structure of the lat¬ 
tice models. We have shown that the self-consistent identity follows from to the Euler- 
Lagrange equation of the pair-correlation functional. This adds a third route for its deriva¬ 
tion to the previous approaches developed by Klamt and Sandler. And it also shows that its 
solution by the Larsen-Rasmussen ansatz corresponds to locating the single minimum of 
the grand potential, while Wilson’s algebraic approach only yields an approximate result. 

Apart from the three investigates models, there exists several more modifications and 
realizations that have not been discussed here. Most often, they differ in their approach 
to solve the self-consistent equation or to couple further interactions, with the Coulomb 
potential as the most relevant example. Including this long-range interaction, adds the 
Debye-Huckel-Pitzer term to the grand potential, associated with a corresponding change 
of the Euler-Lagrange equation. A consistent implementation of the Coulomb interaction 
therefore has to modify the energy functional as well as the self-consistent equation of the 
lattice. 

The main advantage of the lattice models is their computational efficiency. Solving the 
self-consistent equation of the COSMO-RS model only takes seconds. Whereas the calcu¬ 
lation time in the basis of surface segments is of the order of hours. And the minimization 
of the density functional will take even longer. It is therefore necessary to develope further 
approximations for the hard-particle correlations and the integration of the grand potential. 
As for the COSMO-RS model, it is possible to assume a fixed average particle distance 
for the molecular functional and to reduce the integration over the Euclidean volume to a 
sum over all segment pairings and rotations. 

The possibility to derive the thermodynamic equilibrium from a surface integral illus¬ 
trates the potential advantage of the DET ansatz compared to molecular dynamic and 
Monte Carlo simulations. These two methods apply to the complete phase diagram, allow 
for flexible atomic bonds, and the implementation of boundary conditions. But it is quite 
difficult to reduce this freedom, when one is only interested in a small interval of the phase 
diagram or in specific aspects of the intermolecular properties. Well known examples are 
the solubility of proteins, the contact probability of enzymes, the miscibility of racemic 
mixtures, or the selection of an optimal chiral selector in liquid chromatography. These 
are only some examples where the usage of an optimized density functional might prove 
favorable compared to the explicit ensemble averaging of the free energy. 
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Appendix A. Appendix 

The r-particle correlation functionals are always accompanied by an integration over the 
kinematic measure of translations and rotations. Eor the most common case of 

r = 2, we will now derive an explicit realization using methods from integral geometry 
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EMU. 

Let Zjt denote a n — 1 dimensional, smooth, boundary free, convex, Riemannian mani¬ 
fold, imbedded into M”. Each point G is then related to an orthonormal, 

positively oriented coordinate frame ,..., e®) with the outward pointing normal vec- 

(k.) 

tor ef, , differential basis 6,, and connection forms (Oij 

dp — ; dSi — (0jj6j , dCn — (Onak^a — i ^aj8 •— ^afi ' ('^1) 

Each point of the smooth surface is uniquely related to a tangential plane up to an 

(k) 

axial rotation around €„ . Erom this follows that the tangential planes of two convex 
surfaces Zi, Z 2 , touching in a common point pi = p 2 , also agree up to an axial rotation 
and an inversion of their normal vectors This property remains unchanged 

even when the particles are shifted apart in the normal direction. Surface points of closest 
distance t G M+ and their frames are therefore related by 

P 2 = P\ +ten'’ , en'^ = , ea'^ = for G SO(n - 1) (A2) 

using the index conventions i,j = l,...,n and a,jS = 1,... ,n — 1. 

Having defined the relative coordinate frames for the two particles, we can now write 
the integral 


Jg\ 2 {r\,r 2 )dy\dY 2 = J gi 2 (n-^ 2 )^ 712^7 = Evol(SO(n)) J “^ 2 )fi? 7 i 2 (A3) 


in a comoving ^712 and a reference system dy. The integral of the latter can be carried out, 
contributing the volume V = vol(M") and the volume of the group SO(n). The analogous 
transformation in the representation of base forms corresponds to the shift 


dyidy2= A A^/' 

k=l,2 i 


(k) 


A "A = A "A A A 


l(l) fl(2) 


(1) A o(2) 


,(2) 


i<j 


i<i 


i<j 


as can be seen by expanding the skew-symmetric product and setting dy = dy 2 for the 
reference system. 

The trivial contribution dy will be ignored in the following, leaving us with the trans¬ 
formation of dyi 2 - To simplify the calculation, observe that the translation of Z 2 can also 

be written as p 2 + t2^A = P 1 + L A' Aor any 1 1 , t 2 G and ti+t 2 = t. This allows to first 
determine WeyTs half-tube surface Z(t) : p{t) = p + tCn for Zi, Z 2 separately and then to 

derive their Minkowski sum Zi(ti) 0Z2(f2) atpi(ri) = P2{t2) EUSIl- 

In the first step, we determine the differential forms di{t),COij{t) of the half-tube Z(t) by 
differentiating each point p{t) = p + te„ 

dp{t) = dp + indt+t(Ona (A5) 


and separating their components into the directions of and ia 

6n{t) = 6„+dt, Oa{t) = ^a+t(0,ia +th(^p)6p . (A 6 ) 

The new basis also determines the connection forms, as the orthonormal vectors ea and 
their differentials are invariant under translations 


(Ona{t) — tttna 




(A7) 
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Inserting (A 61 , finally yields the curvature matrix and its inverse for the half-tuhe £(?) 

^a /3 (0 ~ ^0:7 (^/jS A t/lyjs) , (0 ~ ^ ^a /3 “h(-^8) 


The second step requires to determine the differential volume element for the domain 
Zi(fi) ©£ 2 ( 12)5 covered hy Zi(li) while circling £ 2 ( 12 )- To write <i 7 i 2 in a common 
coordinate frame we use the transformation ( |A2| ), which relates the connection forms 
cojia = of the two particles at their intersection point pi{ti) = piih) and also 

defines the transformation of their basis forms. Using £ 2 ( 12 ) as reference system, the 
forms of £1 (ti) are rotated into the new coordinate frame 


fl(2) _ T (2) .,(2) 


;(!)„ „( 2 ) 


(A9) 


Inserting this result into (A4), together with the transformation of the normal compo¬ 
nent = —dt and the Jacohi determinant 7 = — 1, yields the reduced kinematic 

measure 


' i<j 

= f\{uX^^\t,)u-^ +X^'^\t2))apCof^ r^dt f\coi^^ A ^al 

a a a<p 

= det(A^^^ +t5 + u^^X^'^'^u) ^da\ Adcj2AdSO{n — 1) Adt , 


where we introduced the unit matrix d, the Gaussian curvature AaCOna = K(}do, the dif¬ 
ferential surface element do, used the orthonormal property det (u) = 1 and (A 81 to write 
the final result in a symmetric form. 
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